#ifndef ORTHAN_HPP
#define ORTHAN_HPP

#include "Density.hpp"
#include <armadillo>
#include <iostream>
#include "utils.hpp"

using namespace arma;
using namespace std;

namespace Orthant{




class Orthant
{
	public:
		Orthant(mat A, mat g,int M){ 
			_A=A;
			_g=g;
			_p=A.n_rows;
			d= new boost::math::normal_distribution<>(0,1); 
			_gamman=0;
			vect.resize(_p,M);
  		}
		~Orthant()
		{
			delete d;
		}
		inline double Likelihood(mat theta1)
		{
			double sum=0;
			for(int i=0;i<_p;i++)
			{
				sum+=log(pdf(*d,theta1(i,0)));

			}
			return sum;
		}

		double Const(mat theta, double gamma, int Bool,int j=1)
		{
			if(!Bool)
			{
			 	vect(span::all,j)=_A*theta-_g.t();	
			}
			double violation=0;
			for(int i=0;i<_p;i++)
			{
				if(vect(i,j)<log(gamma))
				{
					violation+=1;					
				}
			}	
			if(violation>0)
			{
				return -numeric_limits<double>::infinity();
			}else{
				return 0;
			}
	
		}
		double Density_utn(mat theta)
		{
			return Likelihood(theta)+Const(theta,_gamman,0);
		}
		void Weight(mat theta,double *W)
		{
			int M=theta.n_cols;
			for(int i=0;i<M;i++)
			{
				if(W[i]==0)
				{
					if(Const(theta.col(i),_gamman,1,i)==0)
					{
						W[i]=0;
					}else{
						W[i]=-numeric_limits<double>::infinity();
					}
				}
				//W[i]=_Temp(i,0);
				

			}
		}
		double* Eta_star(double gamma, double *W, mat theta, int Bool)
		{
			int M=theta.n_cols;
			_Temp.resize(M,1);
			double *foo=new double[M];
			gamma=_gamman+gamma;
			for(int i=0;i<M;i++)
			{
				if(W[i]==0)
				{
					if(Const(theta.col(i),gamma,Bool,i)==0)
					{
						foo[i]=0;
					}else{
						foo[i]=-numeric_limits<double>::infinity();
					}
				}else{
					foo[i]=W[i];
				}
				_Temp(i,0)=foo[i];
	
			}

			return foo;

		}


		inline void Set_gamma(double gamman){
			_gamma_n1=_gamman;
			_gamman=gamman;
			
		}
		inline double Get_gamma(void){return _gamman;}
		inline int Get_p(void){return _p;}
		boost::math::normal_distribution<>* Get_d(void){return d;}
	private:
		mat _A;
		mat _g;
		mat vect;
		double _gamman;
		double _gamma_n1;
		int _p;
		mat _Temp;
		boost::math::normal_distribution<> *d;

};

} // end namespace
#endif
